– ACTAN.FOR –
FUNCTION ACTAN(SINX,COSX)
COMMON/C2/DE2RA,PI,PIO2,TWOPI,X3PIO2
ACTAN=0.
IF (COSX.EQ.0. ) GO TO 5
IF (COSX.GT.0. ) GO TO 1
ACTAN=PI
GO TO 7
1 IF (SINX.EQ.0. ) GO TO 8
IF (SINX.GT.0. ) GO TO 7
ACTAN=TWOPI
GO TO 7
5 IF (SINX.EQ.0. ) GO TO 8
IF (SINX.GT.0. ) GO TO 6
ACTAN=X3PIO2
GO TO 8
6 ACTAN=PIO2
GO TO 8
7 TEMP=SINX/COSX
ACTAN=ACTAN+ATAN(TEMP)
8 RETURN
END
– DEEP.FOR –
* DEEP SPACE 31 OCT 80
SUBROUTINE DEEP
COMMON/E1/XMO,XNODEO,OMEGAO,EO,XINCL,XNO,XNDT2O,
1 XNDD6O,BSTAR,X,Y,Z,XDOT,YDOT,ZDOT,EPOCH,DS50
COMMON/C1/CK2,CK4,E6A,QOMS2T,S,TOTHRD,
1 XJ3,XKE,XKMPER,XMNPDA,AE
COMMON/C2/DE2RA,PI,PIO2,TWOPI,X3PIO2
DOUBLE PRECISION EPOCH, DS50
DOUBLE PRECISION
* DAY,PREEP,XNODCE,ATIME,DELT,SAVTSN,STEP2,STEPN,STEPP
DATA ZNS, C1SS, ZES/
A 1.19459E-5, 2.9864797E-6, .01675/
DATA ZNL, C1L, ZEL/
A 1.5835218E-4, 4.7968065E-7, .05490/
DATA ZCOSIS, ZSINIS, ZSINGS/
A .91744867, .39785416, -.98088458/
DATA ZCOSGS, ZCOSHS, ZSINHS/
A .1945905, 1.0, 0.0/
DATA Q22,Q31,Q33/1.7891679E-6,2.1460748E-6,2.2123015E-7/
DATA G22,G32/5.7686396,0.95240898/
DATA G44,G52/1.8014998,1.0508330/
DATA G54/4.4108898/
DATA ROOT22,ROOT32/1.7891679E-6,3.7393792E-7/
DATA ROOT44,ROOT52/7.3636953E-9,1.1428639E-7/
DATA ROOT54/2.1765803E-9/
DATA THDT/4.3752691E-3/
* ENTRANCE FOR DEEP SPACE INITIALIZATION
ENTRY DPINIT(EQSQ,SINIQ,COSIQ,RTEQSQ,AO,COSQ2,SINOMO,COSOMO,
1 BSQ,XLLDOT,OMGDT,XNODOT,XNODP)
THGR=THETAG(EPOCH)
EQ = EO
XNQ = XNODP
AQNV = 1./AO
XQNCL = XINCL
XMAO=XMO
XPIDOT=OMGDT+XNODOT
SINQ = SIN(XNODEO)
COSQ = COS(XNODEO)
OMEGAQ = OMEGAO
* INITIALIZE LUNAR SOLAR TERMS
5 DAY=DS50+18261.5D0
IF (DAY.EQ.PREEP) GO TO 10
PREEP = DAY
XNODCE=4.5236020-9.2422029E-4*DAY
STEM=DSIN (XNODCE)
CTEM=DCOS (XNODCE)
ZCOSIL=.91375164-.03568096*CTEM
ZSINIL=SQRT (1.-ZCOSIL*ZCOSIL)
ZSINHL= .089683511*STEM/ZSINIL
ZCOSHL=SQRT (1.-ZSINHL*ZSINHL)
C=4.7199672+.22997150*DAY
GAM=5.8351514+.0019443680*DAY
ZMOL = FMOD2P(C-GAM)
ZX= .39785416*STEM/ZSINIL
ZY= ZCOSHL*CTEM+0.91744867*ZSINHL*STEM
ZX=ACTAN(ZX,ZY)
ZX=GAM+ZX-XNODCE
ZCOSGL=COS (ZX)
ZSINGL=SIN (ZX)
ZMOS=6.2565837D0+.017201977D0*DAY
ZMOS=FMOD2P(ZMOS)
* DO SOLAR TERMS
10 LS = 0
SAVTSN=1.D20
ZCOSG=ZCOSGS
ZSING=ZSINGS
ZCOSI=ZCOSIS
ZSINI=ZSINIS
ZCOSH=COSQ
ZSINH=SINQ
CC=C1SS
ZN=ZNS
ZE=ZES
ZMO=ZMOS
XNOI=1./XNQ
ASSIGN 30 TO LS
20 A1=ZCOSG*ZCOSH+ZSING*ZCOSI*ZSINH
A3=-ZSING*ZCOSH+ZCOSG*ZCOSI*ZSINH
A7=-ZCOSG*ZSINH+ZSING*ZCOSI*ZCOSH
A8=ZSING*ZSINI
A9=ZSING*ZSINH+ZCOSG*ZCOSI*ZCOSH
A10=ZCOSG*ZSINI
A2= COSIQ*A7+ SINIQ*A8
A4= COSIQ*A9+ SINIQ*A10
A5=- SINIQ*A7+ COSIQ*A8
A6=- SINIQ*A9+ COSIQ*A10
C
X1=A1*COSOMO+A2*SINOMO
X2=A3*COSOMO+A4*SINOMO
X3=-A1*SINOMO+A2*COSOMO
X4=-A3*SINOMO+A4*COSOMO
X5=A5*SINOMO
X6=A6*SINOMO
X7=A5*COSOMO
X8=A6*COSOMO
C
Z31=12.*X1*X1-3.*X3*X3
Z32=24.*X1*X2-6.*X3*X4
Z33=12.*X2*X2-3.*X4*X4
Z1=3.*(A1*A1+A2*A2)+Z31*EQSQ
Z2=6.*(A1*A3+A2*A4)+Z32*EQSQ
Z3=3.*(A3*A3+A4*A4)+Z33*EQSQ
Z11=-6.*A1*A5+EQSQ *(-24.*X1*X7-6.*X3*X5)
Z12=-6.*(A1*A6+A3*A5)+EQSQ *(-24.*(X2*X7+X1*X8)-6.*(X3*X6+X4*X5))
Z13=-6.*A3*A6+EQSQ *(-24.*X2*X8-6.*X4*X6)
Z21=6.*A2*A5+EQSQ *(24.*X1*X5-6.*X3*X7)
Z22=6.*(A4*A5+A2*A6)+EQSQ *(24.*(X2*X5+X1*X6)-6.*(X4*X7+X3*X8))
Z23=6.*A4*A6+EQSQ *(24.*X2*X6-6.*X4*X8)
Z1=Z1+Z1+BSQ*Z31
Z2=Z2+Z2+BSQ*Z32
Z3=Z3+Z3+BSQ*Z33
S3=CC*XNOI
S2=-.5*S3/RTEQSQ
S4=S3*RTEQSQ
S1=-15.*EQ*S4
S5=X1*X3+X2*X4
S6=X2*X3+X1*X4
S7=X2*X4-X1*X3
SE=S1*ZN*S5
SI=S2*ZN*(Z11+Z13)
SL=-ZN*S3*(Z1+Z3-14.-6.*EQSQ)
SGH=S4*ZN*(Z31+Z33-6.)
SH=-ZN*S2*(Z21+Z23)
IF(XQNCL.LT.5.2359877E-2) SH=0.0
EE2=2.*S1*S6
E3=2.*S1*S7
XI2=2.*S2*Z12
XI3=2.*S2*(Z13-Z11)
XL2=-2.*S3*Z2
XL3=-2.*S3*(Z3-Z1)
XL4=-2.*S3*(-21.-9.*EQSQ)*ZE
XGH2=2.*S4*Z32
XGH3=2.*S4*(Z33-Z31)
XGH4=-18.*S4*ZE
XH2=-2.*S2*Z22
XH3=-2.*S2*(Z23-Z21)
GO TO LS
* DO LUNAR TERMS
30 SSE = SE
SSI=SI
SSL=SL
SSH=SH/SINIQ
SSG=SGH-COSIQ*SSH
SE2=EE2
SI2=XI2
SL2=XL2
SGH2=XGH2
SH2=XH2
SE3=E3
SI3=XI3
SL3=XL3
SGH3=XGH3
SH3=XH3
SL4=XL4
SGH4=XGH4
LS=1
ZCOSG=ZCOSGL
ZSING=ZSINGL
ZCOSI=ZCOSIL
ZSINI=ZSINIL
ZCOSH=ZCOSHL*COSQ+ZSINHL*SINQ
ZSINH=SINQ*ZCOSHL-COSQ*ZSINHL
ZN=ZNL
CC=C1L
ZE=ZEL
ZMO=ZMOL
ASSIGN 40 TO LS
GO TO 20
40 SSE = SSE+SE
SSI=SSI+SI
SSL=SSL+SL
SSG=SSG+SGH-COSIQ/SINIQ*SH
SSH=SSH+SH/SINIQ
* GEOPOTENTIAL RESONANCE INITIALIZATION FOR 12 HOUR ORBITS
IRESFL=0
ISYNFL=0
IF(XNQ.LT.(.0052359877).AND.XNQ.GT.(.0034906585)) GO TO 70
IF (XNQ.LT.(8.26E-3) .OR. XNQ.GT.(9.24E-3)) RETURN
IF (EQ.LT.0.5) RETURN
IRESFL =1
EOC=EQ*EQSQ
G201=-.306-(EQ-.64)*.440
IF(EQ.GT.(.65)) GO TO 45
G211=3.616-13.247*EQ+16.290*EQSQ
G310=-19.302+117.390*EQ-228.419*EQSQ+156.591*EOC
G322=-18.9068+109.7927*EQ-214.6334*EQSQ+146.5816*EOC
G410=-41.122+242.694*EQ-471.094*EQSQ+313.953*EOC
G422=-146.407+841.880*EQ-1629.014*EQSQ+1083.435*EOC
G520=-532.114+3017.977*EQ-5740*EQSQ+3708.276*EOC
GO TO 55
45 G211=-72.099+331.819*EQ-508.738*EQSQ+266.724*EOC
G310=-346.844+1582.851*EQ-2415.925*EQSQ+1246.113*EOC
G322=-342.585+1554.908*EQ-2366.899*EQSQ+1215.972*EOC
G410=-1052.797+4758.686*EQ-7193.992*EQSQ+3651.957*EOC
G422=-3581.69+16178.11*EQ-24462.77*EQSQ+12422.52*EOC
IF(EQ.GT.(.715)) GO TO 50
G520=1464.74-4664.75*EQ+3763.64*EQSQ
GO TO 55
50 G520=-5149.66+29936.92*EQ-54087.36*EQSQ+31324.56*EOC
55 IF(EQ.GE.(.7)) GO TO 60
G533=-919.2277+4988.61*EQ-9064.77*EQSQ+5542.21*EOC
G521 = -822.71072+4568.6173*EQ-8491.4146*EQSQ+5337.524*EOC
G532 = -853.666+4690.25*EQ-8624.77*EQSQ+5341.4*EOC
GO TO 65
60 G533=-37995.78+161616.52*EQ-229838.2*EQSQ+109377.94*EOC
G521 = -51752.104+218913.95*EQ-309468.16*EQSQ+146349.42*EOC
G532 = -40023.88+170470.89*EQ-242699.48*EQSQ+115605.82*EOC
65 SINI2=SINIQ*SINIQ
F220=.75*(1.+2.*COSIQ+COSQ2)
F221=1.5*SINI2
F321=1.875*SINIQ*(1.-2.*COSIQ-3.*COSQ2)
F322=-1.875*SINIQ*(1.+2.*COSIQ-3.*COSQ2)
F441=35.*SINI2*F220
F442=39.3750*SINI2*SINI2
F522=9.84375*SINIQ*(SINI2*(1.-2.*COSIQ-5.*COSQ2)
1 +.33333333*(-2.+4.*COSIQ+6.*COSQ2))
F523 = SINIQ*(4.92187512*SINI2*(-2.-4.*COSIQ+10.*COSQ2)
* +6.56250012*(1.+2.*COSIQ-3.*COSQ2))
F542 = 29.53125*SINIQ*(2.-8.*COSIQ+COSQ2*(-12.+8.*COSIQ
* +10.*COSQ2))
F543=29.53125*SINIQ*(-2.-8.*COSIQ+COSQ2*(12.+8.*COSIQ-10.*COSQ2))
XNO2=XNQ*XNQ
AINV2=AQNV*AQNV
TEMP1 = 3.*XNO2*AINV2
TEMP = TEMP1*ROOT22
D2201 = TEMP*F220*G201
D2211 = TEMP*F221*G211
TEMP1 = TEMP1*AQNV
TEMP = TEMP1*ROOT32
D3210 = TEMP*F321*G310
D3222 = TEMP*F322*G322
TEMP1 = TEMP1*AQNV
TEMP = 2.*TEMP1*ROOT44
D4410 = TEMP*F441*G410
D4422 = TEMP*F442*G422
TEMP1 = TEMP1*AQNV
TEMP = TEMP1*ROOT52
D5220 = TEMP*F522*G520
D5232 = TEMP*F523*G532
TEMP = 2.*TEMP1*ROOT54
D5421 = TEMP*F542*G521
D5433 = TEMP*F543*G533
XLAMO = XMAO+XNODEO+XNODEO-THGR-THGR
BFACT = XLLDOT+XNODOT+XNODOT-THDT-THDT
BFACT=BFACT+SSL+SSH+SSH
GO TO 80
* SYNCHRONOUS RESONANCE TERMS INITIALIZATION
70 IRESFL=1
ISYNFL=1
G200=1.0+EQSQ*(-2.5+.8125*EQSQ)
G310=1.0+2.0*EQSQ
G300=1.0+EQSQ*(-6.0+6.60937*EQSQ)
F220=.75*(1.+COSIQ)*(1.+COSIQ)
F311=.9375*SINIQ*SINIQ*(1.+3.*COSIQ)-.75*(1.+COSIQ)
F330=1.+COSIQ
F330=1.875*F330*F330*F330
DEL1=3.*XNQ*XNQ*AQNV*AQNV
DEL2=2.*DEL1*F220*G200*Q22
DEL3=3.*DEL1*F330*G300*Q33*AQNV
DEL1=DEL1*F311*G310*Q31*AQNV
FASX2=.13130908
FASX4=2.8843198
FASX6=.37448087
XLAMO=XMAO+XNODEO+OMEGAO-THGR
BFACT = XLLDOT+XPIDOT-THDT
BFACT=BFACT+SSL+SSG+SSH
80 XFACT=BFACT-XNQ
C
C INITIALIZE INTEGRATOR
C
XLI=XLAMO
XNI=XNQ
ATIME=0.D0
STEPP=720.D0
STEPN=-720.D0
STEP2 = 259200.D0
RETURN
* ENTRANCE FOR DEEP SPACE SECULAR EFFECTS
ENTRY DPSEC(XLL,OMGASM,XNODES,EM,XINC,XN,T)
XLL=XLL+SSL*T
OMGASM=OMGASM+SSG*T
XNODES=XNODES+SSH*T
EM=EO+SSE*T
XINC=XINCL+SSI*T
IF(XINC .GE. 0.) GO TO 90
XINC = -XINC
XNODES = XNODES + PI
OMGASM = OMGASM - PI
90 IF(IRESFL .EQ. 0) RETURN
100 IF (ATIME.EQ.0.D0) GO TO 170
IF(T.GE.(0.D0).AND.ATIME.LT.(0.D0)) GO TO 170
IF(T.LT.(0.D0).AND.ATIME.GE.(0.D0)) GO TO 170
105 IF(DABS(T).GE.DABS(ATIME)) GO TO 120
DELT=STEPP
IF (T.GE.0.D0) DELT = STEPN
110 ASSIGN 100 TO IRET
GO TO 160
120 DELT=STEPN
IF (T.GT.0.D0) DELT = STEPP
125 IF (DABS(T-ATIME).LT.STEPP) GO TO 130
ASSIGN 125 TO IRET
GO TO 160
130 FT = T-ATIME
ASSIGN 140 TO IRETN
GO TO 150
140 XN = XNI+XNDOT*FT+XNDDT*FT*FT*0.5
XL = XLI+XLDOT*FT+XNDOT*FT*FT*0.5
TEMP = -XNODES+THGR+T*THDT
XLL = XL-OMGASM+TEMP
IF (ISYNFL.EQ.0) XLL = XL+TEMP+TEMP
RETURN
C
C DOT TERMS CALCULATED
C
150 IF (ISYNFL.EQ.0) GO TO 152
XNDOT=DEL1*SIN (XLI-FASX2)+DEL2*SIN (2.*(XLI-FASX4))
1 +DEL3*SIN (3.*(XLI-FASX6))
XNDDT = DEL1*COS(XLI-FASX2)
* +2.*DEL2*COS(2.*(XLI-FASX4))
* +3.*DEL3*COS(3.*(XLI-FASX6))
GO TO 154
152 XOMI = OMEGAQ+OMGDT*ATIME
X2OMI = XOMI+XOMI
X2LI = XLI+XLI
XNDOT = D2201*SIN(X2OMI+XLI-G22)
* +D2211*SIN(XLI-G22)
* +D3210*SIN(XOMI+XLI-G32)
* +D3222*SIN(-XOMI+XLI-G32)
* +D4410*SIN(X2OMI+X2LI-G44)
* +D4422*SIN(X2LI-G44)
* +D5220*SIN(XOMI+XLI-G52)
* +D5232*SIN(-XOMI+XLI-G52)
* +D5421*SIN(XOMI+X2LI-G54)
* +D5433*SIN(-XOMI+X2LI-G54)
XNDDT = D2201*COS(X2OMI+XLI-G22)
* +D2211*COS(XLI-G22)
* +D3210*COS(XOMI+XLI-G32)
* +D3222*COS(-XOMI+XLI-G32)
* +D5220*COS(XOMI+XLI-G52)
* +D5232*COS(-XOMI+XLI-G52)
* +2.*(D4410*COS(X2OMI+X2LI-G44)
* +D4422*COS(X2LI-G44)
* +D5421*COS(XOMI+X2LI-G54)
* +D5433*COS(-XOMI+X2LI-G54))
154 XLDOT=XNI+XFACT
XNDDT = XNDDT*XLDOT
GO TO IRETN
C
C INTEGRATOR
C
160 ASSIGN 165 TO IRETN
GO TO 150
165 XLI = XLI+XLDOT*DELT+XNDOT*STEP2
XNI = XNI+XNDOT*DELT+XNDDT*STEP2
ATIME=ATIME+DELT
GO TO IRET
C
C EPOCH RESTART
C
170 IF (T.GE.0.D0) GO TO 175
DELT=STEPN
GO TO 180
175 DELT = STEPP
180 ATIME = 0.D0
XNI=XNQ
XLI=XLAMO
GO TO 125
C
C ENTRANCES FOR LUNAR-SOLAR PERIODICS
C
C
ENTRY DPPER(EM,XINC,OMGASM,XNODES,XLL)
SINIS = SIN(XINC)
COSIS = COS(XINC)
IF (DABS(SAVTSN-T).LT.(30.D0)) GO TO 210
SAVTSN=T
ZM=ZMOS+ZNS*T
205 ZF=ZM+2.*ZES*SIN (ZM)
SINZF=SIN (ZF)
F2=.5*SINZF*SINZF-.25
F3=-.5*SINZF*COS (ZF)
SES=SE2*F2+SE3*F3
SIS=SI2*F2+SI3*F3
SLS=SL2*F2+SL3*F3+SL4*SINZF
SGHS=SGH2*F2+SGH3*F3+SGH4*SINZF
SHS=SH2*F2+SH3*F3
ZM=ZMOL+ZNL*T
ZF=ZM+2.*ZEL*SIN (ZM)
SINZF=SIN (ZF)
F2=.5*SINZF*SINZF-.25
F3=-.5*SINZF*COS (ZF)
SEL=EE2*F2+E3*F3
SIL=XI2*F2+XI3*F3
SLL=XL2*F2+XL3*F3+XL4*SINZF
SGHL=XGH2*F2+XGH3*F3+XGH4*SINZF
SHL=XH2*F2+XH3*F3
PE=SES+SEL
PINC=SIS+SIL
PL=SLS+SLL
210 PGH=SGHS+SGHL
PH=SHS+SHL
XINC = XINC+PINC
EM = EM+PE
IF(XQNCL.LT.(.2)) GO TO 220
GO TO 218
C
C APPLY PERIODICS DIRECTLY
C
218 PH=PH/SINIQ
PGH=PGH-COSIQ*PH
OMGASM=OMGASM+PGH
XNODES=XNODES+PH
XLL = XLL+PL
GO TO 230
C
C APPLY PERIODICS WITH LYDDANE MODIFICATION
C
220 SINOK=SIN(XNODES)
COSOK=COS(XNODES)
ALFDP=SINIS*SINOK
BETDP=SINIS*COSOK
DALF=PH*COSOK+PINC*COSIS*SINOK
DBET=-PH*SINOK+PINC*COSIS*COSOK
ALFDP=ALFDP+DALF
BETDP=BETDP+DBET
XLS = XLL+OMGASM+COSIS*XNODES
DLS=PL+PGH-PINC*XNODES*SINIS
XLS=XLS+DLS
XNODES=ACTAN(ALFDP,BETDP)
XLL = XLL+PL
OMGASM = XLS-XLL-COS(XINC)*XNODES
230 CONTINUE
RETURN
END
– DRIVER.FOR –
* DRIVER 3 NOV 80
* WGS-72 PHYSICAL AND GEOPOTENTIAL CONSTANTS
* CK2= .5*J2*AE**2 CK4=-.375*J4*AE**4
DOUBLE PRECISION EPOCH,DS50
COMMON/E1/XMO,XNODEO,OMEGAO,EO,XINCL,XNO,XNDT2O,XNDD6O,BSTAR,
1 X,Y,Z,XDOT,YDOT,ZDOT,EPOCH,DS50
COMMON/C1/CK2,CK4,E6A,QOMS2T,S,TOTHRD,
1 XJ3,XKE,XKMPER,XMNPDA,AE
COMMON/C2/DE2RA,PI,PIO2,TWOPI,X3PIO2
DATA IHG/1HG/
DATA DE2RA,E6A,PI,PIO2,QO,SO,TOTHRD,TWOPI,X3PIO2,XJ2,XJ3,
1 XJ4,XKE,XKMPER,XMNPDA,AE/.174532925E-1,1.E-6,
2 3.14159265,1.57079633,120.0,78.0,.66666667,
4 6.2831853,4.71238898,1.082616E-3,-.253881E-5,
5 -1.65597E-6,.743669161E-1,6378.135,1440.,1./
DIMENSION ISET(5)
CHARACTER ABUF*80(2)
DATA (ISET(I),I=1,5)/3HSGP,4HSGP4,4HSDP4,4HSGP8,4HSDP8/
* SELECT EPHEMERIS TYPE AND OUTPUT TIMES
CK2=.5*XJ2*AE**2
CK4=-.375*XJ4*AE**4
QOMS2T=((QO-SO)*AE/XKMPER)**4
S=AE*(1.+SO/XKMPER)
2 READ (5,700) IEPT, TS,TF,DELT
IF(IEPT.LE.0) STOP
IDEEP=0
* READ IN MEAN ELEMENTS FROM 2 CARD T(TRANS) OR G(INTERN) FORMAT
READ (5,706) ABUF
DECODE(ABUF(1),707) ITYPE
IF(ITYPE.EQ.IHG) GO TO 5
DECODE (ABUF,702) EPOCH,XNDT2O,XNDD6O,IEXP,BSTAR,IBEXP,XINCL,
1 XNODEO,EO,OMEGAO,XMO,XNO
GO TO 7
5 DECODE(ABUF,701) EPOCH,XMO,XNODEO,OMEGAO,EO,XINCL,XNO,XNDT2O,
1 XNDD6O,IEXP,BSTAR,IBEXP
7 IF(XNO.LE.0.) STOP
WRITE(6,704) ABUF,ISET(IEPT)
IF(IEPT.GT.5) GO TO 900
XNDD6O=XNDD6O*(10.**IEXP)
XNODEO=XNODEO*DE2RA
OMEGAO=OMEGAO*DE2RA
XMO=XMO*DE2RA
XINCL=XINCL*DE2RA
TEMP=TWOPI/XMNPDA/XMNPDA
XNO=XNO*TEMP*XMNPDA
XNDT2O=XNDT2O*TEMP
XNDD6O=XNDD6O*TEMP/XMNPDA
* INPUT CHECK FOR PERIOD VS EPHEMERIS SELECTED
* PERIOD GE 225 MINUTES IS DEEP SPACE
A1=(XKE/XNO)**TOTHRD
TEMP=1.5*CK2*(3.*COS(XINCL)**2-1.)/(1.-EO*EO)**1.5
DEL1=TEMP/(A1*A1)
AO=A1*(1.-DEL1*(.5*TOTHRD+DEL1*(1.+134./81.*DEL1)))
DELO=TEMP/(AO*AO)
XNODP=XNO/(1.+DELO)
IF((TWOPI/XNODP/XMNPDA) .GE. .15625) IDEEP=1
BSTAR=BSTAR*(10.**IBEXP)/AE
TSINCE=TS
IFLAG=1
IF(IDEEP .EQ. 1 .AND. (IEPT .EQ. 1 .OR. IEPT .EQ. 2
1 .OR. IEPT .EQ. 4)) GO TO 800
9 IF(IDEEP .EQ. 0 .AND. (IEPT .EQ. 3 .OR. IEPT .EQ. 5))
1 GO TO 850
10 GO TO (21,22,23,24,25), IEPT
21 CALL SGP(IFLAG,TSINCE)
GO TO 60
22 CALL SGP4(IFLAG,TSINCE)
GO TO 60
23 CALL SDP4(IFLAG,TSINCE)
GO TO 60
24 CALL SGP8(IFLAG,TSINCE)
GO TO 60
25 CALL SDP8(IFLAG,TSINCE)
60 X=X*XKMPER/AE
Y=Y*XKMPER/AE
Z=Z*XKMPER/AE
XDOT=XDOT*XKMPER/AE*XMNPDA/86400.
YDOT=YDOT*XKMPER/AE*XMNPDA/86400.
ZDOT=ZDOT*XKMPER/AE*XMNPDA/86400.
WRITE(6,705) TSINCE,X,Y,Z,XDOT,YDOT,ZDOT
TSINCE=TSINCE+DELT
IF(ABS(TSINCE) .GT. ABS(TF)) GO TO 2
GO TO 10
700 FORMAT(I1,3F10.0)
701 FORMAT(29X,D14.8,1X,3F8.4,/,6X,F8.7,F8.4,1X,2F11.9,1X,F6.5,I2,
1 4X,F8.7,I2)
702 FORMAT(18X,D14.8,1X,F10.8,2(1X,F6.5,I2),/,7X,2(1X,F8.4),1X,
1 F7.7,2(1X,F8.4),1X,F11.8)
703 FORMAT(79X,A1)
704 FORMAT(1H1,A80,/,1X,A80,//,1X,A4,7H TSINCE,
1 14X,1HX,16X,1HY,16X,1HZ,14X,
1 4HXDOT,13X,4HYDOT,13X,4HZDOT,//)
705 FORMAT(7F17.8)
706 FORMAT(A80)
707 FORMAT(79X,A1)
930 FORMAT("SHOULD USE DEEP SPACE EPHEMERIS")
940 FORMAT("SHOULD USE NEAR EARTH EPHEMERIS")
950 FORMAT("EPHEMERIS NUMBER",I2," NOT LEGAL, WILL SKIP THIS CASE")
800 WRITE(6,930)
GO TO 9
850 WRITE(6,940)
GO TO 10
900 WRITE(6,950) IEPT
GO TO 2
END
– FMOD2P.FOR –
FUNCTION FMOD2P(X)
COMMON/C2/DE2RA,PI,PIO2,TWOPI,X3PIO2
FMOD2P=X
I=FMOD2P/TWOPI
FMOD2P=FMOD2P-I*TWOPI
IF(FMOD2P.LT.0) FMOD2P=FMOD2P+TWOPI
RETURN
END
– SGP.FOR –
* SGP 31 OCT 80
SUBROUTINE SGP(IFLAG,TSINCE)
COMMON/E1/XMO,XNODEO,OMEGAO,EO,XINCL,XNO,XNDT2O,XNDD6O,BSTAR,
1 X,Y,Z,XDOT,YDOT,ZDOT,EPOCH,DS50
COMMON/C1/CK2,CK4,E6A,QOMS2T,S,TOTHRD,
1 XJ3,XKE,XKMPER,XMNPDA,AE
DOUBLE PRECISION EPOCH, DS50
IF(IFLAG.EQ.0) GO TO 19
* INITIALIZATION
C1= CK2*1.5
C2= CK2/4.0
C3= CK2/2.0
C4= XJ3*AE**3/(4.0*CK2)
COSIO=COS(XINCL)
SINIO=SIN(XINCL)
A1=(XKE/XNO)**TOTHRD
D1= C1/A1/A1*(3.*COSIO*COSIO-1.)/(1.-EO*EO)**1.5
AO=A1*(1.-1./3.*D1-D1*D1-134./81.*D1*D1*D1)
PO=AO*(1.-EO*EO)
QO=AO*(1.-EO)
XLO=XMO+OMEGAO+XNODEO
D1O= C3 *SINIO*SINIO
D2O= C2 *(7.*COSIO*COSIO-1.)
D3O=C1*COSIO
D4O=D3O*SINIO
PO2NO=XNO/(PO*PO)
OMGDT=C1*PO2NO*(5.*COSIO*COSIO-1.)
XNODOT=-2.*D3O*PO2NO
C5=.5*C4*SINIO*(3.+5.*COSIO)/(1.+COSIO)
C6=C4*SINIO
IFLAG=0
* UPDATE FOR SECULAR GRAVITY AND ATMOSPHERIC DRAG
19 A=XNO+(2.*XNDT2O+3.*XNDD6O*TSINCE)*TSINCE
A=AO*(XNO/A)**TOTHRD
E=E6A
IF(A.GT.QO) E=1.-QO/A
P=A*(1.-E*E)
XNODES= XNODEO+XNODOT*TSINCE
OMGAS= OMEGAO+OMGDT*TSINCE
XLS=FMOD2P(XLO+(XNO+OMGDT+XNODOT+(XNDT2O+XNDD6O*TSINCE)*
1 TSINCE)*TSINCE)
* LONG PERIOD PERIODICS
AXNSL=E*COS(OMGAS)
AYNSL=E*SIN(OMGAS)-C6/P
XL=FMOD2P(XLS-C5/P*AXNSL)
* SOLVE KEPLERS EQUATION
U=FMOD2P(XL-XNODES)
ITEM3=0
EO1=U
TEM5=1.
20 SINEO1=SIN(EO1)
COSEO1=COS(EO1)
IF(ABS(TEM5).LT.E6A) GO TO 30
IF(ITEM3.GE.10) GO TO 30
ITEM3=ITEM3+1
TEM5=1.-COSEO1*AXNSL-SINEO1*AYNSL
TEM5=(U-AYNSL*COSEO1+AXNSL*SINEO1-EO1)/TEM5
TEM2=ABS(TEM5)
IF(TEM2.GT.1.) TEM5=TEM2/TEM5
EO1=EO1+TEM5
GO TO 20
* SHORT PERIOD PRELIMINARY QUANTITIES
30 ECOSE=AXNSL*COSEO1+AYNSL*SINEO1
ESINE=AXNSL*SINEO1-AYNSL*COSEO1
EL2=AXNSL*AXNSL+AYNSL*AYNSL
PL=A*(1.-EL2)
PL2=PL*PL
R=A*(1.-ECOSE)
RDOT=XKE*SQRT(A)/R*ESINE
RVDOT=XKE*SQRT(PL)/R
TEMP=ESINE/(1.+SQRT(1.-EL2))
SINU=A/R*(SINEO1-AYNSL-AXNSL*TEMP)
COSU=A/R*(COSEO1-AXNSL+AYNSL*TEMP)
SU=ACTAN(SINU,COSU)
* UPDATE FOR SHORT PERIODICS
SIN2U=(COSU+COSU)*SINU
COS2U=1.-2.*SINU*SINU
RK=R+D1O/PL*COS2U
UK=SU-D2O/PL2*SIN2U
XNODEK=XNODES+D3O*SIN2U/PL2
XINCK =XINCL+D4O/PL2*COS2U
* ORIENTATION VECTORS
SINUK=SIN(UK)
COSUK=COS(UK)
SINNOK=SIN(XNODEK)
COSNOK=COS(XNODEK)
SINIK=SIN(XINCK)
COSIK=COS(XINCK)
XMX=-SINNOK*COSIK
XMY=COSNOK*COSIK
UX=XMX*SINUK+COSNOK*COSUK
UY=XMY*SINUK+SINNOK*COSUK
UZ=SINIK*SINUK
VX=XMX*COSUK-COSNOK*SINUK
VY=XMY*COSUK-SINNOK*SINUK
VZ=SINIK*COSUK
* POSITION AND VELOCITY
X=RK*UX
Y=RK*UY
Z=RK*UZ
XDOT=RDOT*UX
YDOT=RDOT*UY
ZDOT=RDOT*UZ
XDOT=RVDOT*VX+XDOT
YDOT=RVDOT*VY+YDOT
ZDOT=RVDOT*VZ+ZDOT
RETURN
END
– SGP4.FOR –
* SGP4 3 NOV 80
SUBROUTINE SGP4(IFLAG,TSINCE)
COMMON/E1/XMO,XNODEO,OMEGAO,EO,XINCL,XNO,XNDT2O,
1 XNDD6O,BSTAR,X,Y,Z,XDOT,YDOT,ZDOT,EPOCH,DS50
COMMON/C1/CK2,CK4,E6A,QOMS2T,S,TOTHRD,
1 XJ3,XKE,XKMPER,XMNPDA,AE
DOUBLE PRECISION EPOCH, DS50
IF (IFLAG .EQ. 0) GO TO 100
* RECOVER ORIGINAL MEAN MOTION (XNODP) AND SEMIMAJOR AXIS (AODP)
* FROM INPUT ELEMENTS
A1=(XKE/XNO)**TOTHRD
COSIO=COS(XINCL)
THETA2=COSIO*COSIO
X3THM1=3.*THETA2-1.
EOSQ=EO*EO
BETAO2=1.-EOSQ
BETAO=SQRT(BETAO2)
DEL1=1.5*CK2*X3THM1/(A1*A1*BETAO*BETAO2)
AO=A1*(1.-DEL1*(.5*TOTHRD+DEL1*(1.+134./81.*DEL1)))
DELO=1.5*CK2*X3THM1/(AO*AO*BETAO*BETAO2)
XNODP=XNO/(1.+DELO)
AODP=AO/(1.-DELO)
* INITIALIZATION
* FOR PERIGEE LESS THAN 220 KILOMETERS, THE ISIMP FLAG IS SET AND
* THE EQUATIONS ARE TRUNCATED TO LINEAR VARIATION IN SQRT A AND
* QUADRATIC VARIATION IN MEAN ANOMALY. ALSO, THE C3 TERM, THE
* DELTA OMEGA TERM, AND THE DELTA M TERM ARE DROPPED.
ISIMP=0
IF((AODP*(1.-EO)/AE) .LT. (220./XKMPER+AE)) ISIMP=1
* FOR PERIGEE BELOW 156 KM, THE VALUES OF
* S AND QOMS2T ARE ALTERED
S4=S
QOMS24=QOMS2T
PERIGE=(AODP*(1.-EO)-AE)*XKMPER
IF(PERIGE .GE. 156.) GO TO 10
S4=PERIGE-78.
IF(PERIGE .GT. 98.) GO TO 9
S4=20.
9 QOMS24=((120.-S4)*AE/XKMPER)**4
S4=S4/XKMPER+AE
10 PINVSQ=1./(AODP*AODP*BETAO2*BETAO2)
TSI=1./(AODP-S4)
ETA=AODP*EO*TSI
ETASQ=ETA*ETA
EETA=EO*ETA
PSISQ=ABS(1.-ETASQ)
COEF=QOMS24*TSI**4
COEF1=COEF/PSISQ**3.5
C2=COEF1*XNODP*(AODP*(1.+1.5*ETASQ+EETA*(4.+ETASQ))+.75*
1 CK2*TSI/PSISQ*X3THM1*(8.+3.*ETASQ*(8.+ETASQ)))
C1=BSTAR*C2
SINIO=SIN(XINCL)
A3OVK2=-XJ3/CK2*AE**3
C3=COEF*TSI*A3OVK2*XNODP*AE*SINIO/EO
X1MTH2=1.-THETA2
C4=2.*XNODP*COEF1*AODP*BETAO2*(ETA*
1 (2.+.5*ETASQ)+EO*(.5+2.*ETASQ)-2.*CK2*TSI/
2 (AODP*PSISQ)*(-3.*X3THM1*(1.-2.*EETA+ETASQ*
3 (1.5-.5*EETA))+.75*X1MTH2*(2.*ETASQ-EETA*
4 (1.+ETASQ))*COS(2.*OMEGAO)))
C5=2.*COEF1*AODP*BETAO2*(1.+2.75*(ETASQ+EETA)+EETA*ETASQ)
THETA4=THETA2*THETA2
TEMP1=3.*CK2*PINVSQ*XNODP
TEMP2=TEMP1*CK2*PINVSQ
TEMP3=1.25*CK4*PINVSQ*PINVSQ*XNODP
XMDOT=XNODP+.5*TEMP1*BETAO*X3THM1+.0625*TEMP2*BETAO*
1 (13.-78.*THETA2+137.*THETA4)
X1M5TH=1.-5.*THETA2
OMGDOT=-.5*TEMP1*X1M5TH+.0625*TEMP2*(7.-114.*THETA2+
1 395.*THETA4)+TEMP3*(3.-36.*THETA2+49.*THETA4)
XHDOT1=-TEMP1*COSIO
XNODOT=XHDOT1+(.5*TEMP2*(4.-19.*THETA2)+2.*TEMP3*(3.-
1 7.*THETA2))*COSIO
OMGCOF=BSTAR*C3*COS(OMEGAO)
XMCOF=-TOTHRD*COEF*BSTAR*AE/EETA
XNODCF=3.5*BETAO2*XHDOT1*C1
T2COF=1.5*C1
XLCOF=.125*A3OVK2*SINIO*(3.+5.*COSIO)/(1.+COSIO)
AYCOF=.25*A3OVK2*SINIO
DELMO=(1.+ETA*COS(XMO))**3
SINMO=SIN(XMO)
X7THM1=7.*THETA2-1.
IF(ISIMP .EQ. 1) GO TO 90
C1SQ=C1*C1
D2=4.*AODP*TSI*C1SQ
TEMP=D2*TSI*C1/3.
D3=(17.*AODP+S4)*TEMP
D4=.5*TEMP*AODP*TSI*(221.*AODP+31.*S4)*C1
T3COF=D2+2.*C1SQ
T4COF=.25*(3.*D3+C1*(12.*D2+10.*C1SQ))
T5COF=.2*(3.*D4+12.*C1*D3+6.*D2*D2+15.*C1SQ*(
1 2.*D2+C1SQ))
90 IFLAG=0
* UPDATE FOR SECULAR GRAVITY AND ATMOSPHERIC DRAG
100 XMDF=XMO+XMDOT*TSINCE
OMGADF=OMEGAO+OMGDOT*TSINCE
XNODDF=XNODEO+XNODOT*TSINCE
OMEGA=OMGADF
XMP=XMDF
TSQ=TSINCE*TSINCE
XNODE=XNODDF+XNODCF*TSQ
TEMPA=1.-C1*TSINCE
TEMPE=BSTAR*C4*TSINCE
TEMPL=T2COF*TSQ
IF(ISIMP .EQ. 1) GO TO 110
DELOMG=OMGCOF*TSINCE
DELM=XMCOF*((1.+ETA*COS(XMDF))**3-DELMO)
TEMP=DELOMG+DELM
XMP=XMDF+TEMP
OMEGA=OMGADF-TEMP
TCUBE=TSQ*TSINCE
TFOUR=TSINCE*TCUBE
TEMPA=TEMPA-D2*TSQ-D3*TCUBE-D4*TFOUR
TEMPE=TEMPE+BSTAR*C5*(SIN(XMP)-SINMO)
TEMPL=TEMPL+T3COF*TCUBE+
1 TFOUR*(T4COF+TSINCE*T5COF)
110 A=AODP*TEMPA**2
E=EO-TEMPE
XL=XMP+OMEGA+XNODE+XNODP*TEMPL
BETA=SQRT(1.-E*E)
XN=XKE/A**1.5
* LONG PERIOD PERIODICS
AXN=E*COS(OMEGA)
TEMP=1./(A*BETA*BETA)
XLL=TEMP*XLCOF*AXN
AYNL=TEMP*AYCOF
XLT=XL+XLL
AYN=E*SIN(OMEGA)+AYNL
* SOLVE KEPLERS EQUATION
CAPU=FMOD2P(XLT-XNODE)
TEMP2=CAPU
DO 130 I=1,10
SINEPW=SIN(TEMP2)
COSEPW=COS(TEMP2)
TEMP3=AXN*SINEPW
TEMP4=AYN*COSEPW
TEMP5=AXN*COSEPW
TEMP6=AYN*SINEPW
EPW=(CAPU-TEMP4+TEMP3-TEMP2)/(1.-TEMP5-TEMP6)+TEMP2
IF(ABS(EPW-TEMP2) .LE. E6A) GO TO 140
130 TEMP2=EPW
* SHORT PERIOD PRELIMINARY QUANTITIES
140 ECOSE=TEMP5+TEMP6
ESINE=TEMP3-TEMP4
ELSQ=AXN*AXN+AYN*AYN
TEMP=1.-ELSQ
PL=A*TEMP
R=A*(1.-ECOSE)
TEMP1=1./R
RDOT=XKE*SQRT(A)*ESINE*TEMP1
RFDOT=XKE*SQRT(PL)*TEMP1
TEMP2=A*TEMP1
BETAL=SQRT(TEMP)
TEMP3=1./(1.+BETAL)
COSU=TEMP2*(COSEPW-AXN+AYN*ESINE*TEMP3)
SINU=TEMP2*(SINEPW-AYN-AXN*ESINE*TEMP3)
U=ACTAN(SINU,COSU)
SIN2U=2.*SINU*COSU
COS2U=2.*COSU*COSU-1.
TEMP=1./PL
TEMP1=CK2*TEMP
TEMP2=TEMP1*TEMP
* UPDATE FOR SHORT PERIODICS
RK=R*(1.-1.5*TEMP2*BETAL*X3THM1)+.5*TEMP1*X1MTH2*COS2U
UK=U-.25*TEMP2*X7THM1*SIN2U
XNODEK=XNODE+1.5*TEMP2*COSIO*SIN2U
XINCK=XINCL+1.5*TEMP2*COSIO*SINIO*COS2U
RDOTK=RDOT-XN*TEMP1*X1MTH2*SIN2U
RFDOTK=RFDOT+XN*TEMP1*(X1MTH2*COS2U+1.5*X3THM1)
* ORIENTATION VECTORS
SINUK=SIN(UK)
COSUK=COS(UK)
SINIK=SIN(XINCK)
COSIK=COS(XINCK)
SINNOK=SIN(XNODEK)
COSNOK=COS(XNODEK)
XMX=-SINNOK*COSIK
XMY=COSNOK*COSIK
UX=XMX*SINUK+COSNOK*COSUK
UY=XMY*SINUK+SINNOK*COSUK
UZ=SINIK*SINUK
VX=XMX*COSUK-COSNOK*SINUK
VY=XMY*COSUK-SINNOK*SINUK
VZ=SINIK*COSUK
* POSITION AND VELOCITY
X=RK*UX
Y=RK*UY
Z=RK*UZ
XDOT=RDOTK*UX+RFDOTK*VX
YDOT=RDOTK*UY+RFDOTK*VY
ZDOT=RDOTK*UZ+RFDOTK*VZ
RETURN
END
– SDP4.FOR –
* SDP4 3 NOV 80
SUBROUTINE SDP4(IFLAG,TSINCE)
COMMON/E1/XMO,XNODEO,OMEGAO,EO,XINCL,XNO,XNDT2O,
1 XNDD6O,BSTAR,X,Y,Z,XDOT,YDOT,ZDOT,EPOCH,DS50
COMMON/C1/CK2,CK4,E6A,QOMS2T,S,TOTHRD,
1 XJ3,XKE,XKMPER,XMNPDA,AE
DOUBLE PRECISION EPOCH, DS50
IF (IFLAG .EQ. 0) GO TO 100
* RECOVER ORIGINAL MEAN MOTION (XNODP) AND SEMIMAJOR AXIS (AODP)
* FROM INPUT ELEMENTS
A1=(XKE/XNO)**TOTHRD
COSIO=COS(XINCL)
THETA2=COSIO*COSIO
X3THM1=3.*THETA2-1.
EOSQ=EO*EO
BETAO2=1.-EOSQ
BETAO=SQRT(BETAO2)
DEL1=1.5*CK2*X3THM1/(A1*A1*BETAO*BETAO2)
AO=A1*(1.-DEL1*(.5*TOTHRD+DEL1*(1.+134./81.*DEL1)))
DELO=1.5*CK2*X3THM1/(AO*AO*BETAO*BETAO2)
XNODP=XNO/(1.+DELO)
AODP=AO/(1.-DELO)
* INITIALIZATION
* FOR PERIGEE BELOW 156 KM, THE VALUES OF
* S AND QOMS2T ARE ALTERED
S4=S
QOMS24=QOMS2T
PERIGE=(AODP*(1.-EO)-AE)*XKMPER
IF(PERIGE .GE. 156.) GO TO 10
S4=PERIGE-78.
IF(PERIGE .GT. 98.) GO TO 9
S4=20.
9 QOMS24=((120.-S4)*AE/XKMPER)**4
S4=S4/XKMPER+AE
10 PINVSQ=1./(AODP*AODP*BETAO2*BETAO2)
SING=SIN(OMEGAO)
COSG=COS(OMEGAO)
TSI=1./(AODP-S4)
ETA=AODP*EO*TSI
ETASQ=ETA*ETA
EETA=EO*ETA
PSISQ=ABS(1.-ETASQ)
COEF=QOMS24*TSI**4
COEF1=COEF/PSISQ**3.5
C2=COEF1*XNODP*(AODP*(1.+1.5*ETASQ+EETA*(4.+ETASQ))+.75*
1 CK2*TSI/PSISQ*X3THM1*(8.+3.*ETASQ*(8.+ETASQ)))
C1=BSTAR*C2
SINIO=SIN(XINCL)
A3OVK2=-XJ3/CK2*AE**3
X1MTH2=1.-THETA2
C4=2.*XNODP*COEF1*AODP*BETAO2*(ETA*
1 (2.+.5*ETASQ)+EO*(.5+2.*ETASQ)-2.*CK2*TSI/
2 (AODP*PSISQ)*(-3.*X3THM1*(1.-2.*EETA+ETASQ*
3 (1.5-.5*EETA))+.75*X1MTH2*(2.*ETASQ-EETA*
4 (1.+ETASQ))*COS(2.*OMEGAO)))
THETA4=THETA2*THETA2
TEMP1=3.*CK2*PINVSQ*XNODP
TEMP2=TEMP1*CK2*PINVSQ
TEMP3=1.25*CK4*PINVSQ*PINVSQ*XNODP
XMDOT=XNODP+.5*TEMP1*BETAO*X3THM1+.0625*TEMP2*BETAO*
1 (13.-78.*THETA2+137.*THETA4)
X1M5TH=1.-5.*THETA2
OMGDOT=-.5*TEMP1*X1M5TH+.0625*TEMP2*(7.-114.*THETA2+
1 395.*THETA4)+TEMP3*(3.-36.*THETA2+49.*THETA4)
XHDOT1=-TEMP1*COSIO
XNODOT=XHDOT1+(.5*TEMP2*(4.-19.*THETA2)+2.*TEMP3*(3.-
1 7.*THETA2))*COSIO
XNODCF=3.5*BETAO2*XHDOT1*C1
T2COF=1.5*C1
XLCOF=.125*A3OVK2*SINIO*(3.+5.*COSIO)/(1.+COSIO)
AYCOF=.25*A3OVK2*SINIO
X7THM1=7.*THETA2-1.
90 IFLAG=0
CALL DPINIT(EOSQ,SINIO,COSIO,BETAO,AODP,THETA2,
1 SING,COSG,BETAO2,XMDOT,OMGDOT,XNODOT,XNODP)
* UPDATE FOR SECULAR GRAVITY AND ATMOSPHERIC DRAG
100 XMDF=XMO+XMDOT*TSINCE
OMGADF=OMEGAO+OMGDOT*TSINCE
XNODDF=XNODEO+XNODOT*TSINCE
TSQ=TSINCE*TSINCE
XNODE=XNODDF+XNODCF*TSQ
TEMPA=1.-C1*TSINCE
TEMPE=BSTAR*C4*TSINCE
TEMPL=T2COF*TSQ
XN=XNODP
CALL DPSEC(XMDF,OMGADF,XNODE,EM,XINC,XN,TSINCE)
A=(XKE/XN)**TOTHRD*TEMPA**2
E=EM-TEMPE
XMAM=XMDF+XNODP*TEMPL
CALL DPPER(E,XINC,OMGADF,XNODE,XMAM)
XL=XMAM+OMGADF+XNODE
BETA=SQRT(1.-E*E)
XN=XKE/A**1.5
* LONG PERIOD PERIODICS
AXN=E*COS(OMGADF)
TEMP=1./(A*BETA*BETA)
XLL=TEMP*XLCOF*AXN
AYNL=TEMP*AYCOF
XLT=XL+XLL
AYN=E*SIN(OMGADF)+AYNL
* SOLVE KEPLERS EQUATION
CAPU=FMOD2P(XLT-XNODE)
TEMP2=CAPU
DO 130 I=1,10
SINEPW=SIN(TEMP2)
COSEPW=COS(TEMP2)
TEMP3=AXN*SINEPW
TEMP4=AYN*COSEPW
TEMP5=AXN*COSEPW
TEMP6=AYN*SINEPW
EPW=(CAPU-TEMP4+TEMP3-TEMP2)/(1.-TEMP5-TEMP6)+TEMP2
IF(ABS(EPW-TEMP2) .LE. E6A) GO TO 140
130 TEMP2=EPW
* SHORT PERIOD PRELIMINARY QUANTITIES
140 ECOSE=TEMP5+TEMP6
ESINE=TEMP3-TEMP4
ELSQ=AXN*AXN+AYN*AYN
TEMP=1.-ELSQ
PL=A*TEMP
R=A*(1.-ECOSE)
TEMP1=1./R
RDOT=XKE*SQRT(A)*ESINE*TEMP1
RFDOT=XKE*SQRT(PL)*TEMP1
TEMP2=A*TEMP1
BETAL=SQRT(TEMP)
TEMP3=1./(1.+BETAL)
COSU=TEMP2*(COSEPW-AXN+AYN*ESINE*TEMP3)
SINU=TEMP2*(SINEPW-AYN-AXN*ESINE*TEMP3)
U=ACTAN(SINU,COSU)
SIN2U=2.*SINU*COSU
COS2U=2.*COSU*COSU-1.
TEMP=1./PL
TEMP1=CK2*TEMP
TEMP2=TEMP1*TEMP
* UPDATE FOR SHORT PERIODICS
RK=R*(1.-1.5*TEMP2*BETAL*X3THM1)+.5*TEMP1*X1MTH2*COS2U
UK=U-.25*TEMP2*X7THM1*SIN2U
XNODEK=XNODE+1.5*TEMP2*COSIO*SIN2U
XINCK=XINC+1.5*TEMP2*COSIO*SINIO*COS2U
RDOTK=RDOT-XN*TEMP1*X1MTH2*SIN2U
RFDOTK=RFDOT+XN*TEMP1*(X1MTH2*COS2U+1.5*X3THM1)
* ORIENTATION VECTORS
SINUK=SIN(UK)
COSUK=COS(UK)
SINIK=SIN(XINCK)
COSIK=COS(XINCK)
SINNOK=SIN(XNODEK)
COSNOK=COS(XNODEK)
XMX=-SINNOK*COSIK
XMY=COSNOK*COSIK
UX=XMX*SINUK+COSNOK*COSUK
UY=XMY*SINUK+SINNOK*COSUK
UZ=SINIK*SINUK
VX=XMX*COSUK-COSNOK*SINUK
VY=XMY*COSUK-SINNOK*SINUK
VZ=SINIK*COSUK
* POSITION AND VELOCITY
X=RK*UX
Y=RK*UY
Z=RK*UZ
XDOT=RDOTK*UX+RFDOTK*VX
YDOT=RDOTK*UY+RFDOTK*VY
ZDOT=RDOTK*UZ+RFDOTK*VZ
RETURN
END
– SGP8.FOR –
* SGP8 14 NOV 80
SUBROUTINE SGP8(IFLAG,TSINCE)
COMMON/E1/XMO,XNODEO,OMEGAO,EO,XINCL,XNO,XNDT2O,
1 XNDD6O,BSTAR,X,Y,Z,XDOT,YDOT,ZDOT,EPOCH,DS50
COMMON/C1/CK2,CK4,E6A,QOMS2T,S,TOTHRD,
1 XJ3,XKE,XKMPER,XMNPDA,AE
DOUBLE PRECISION EPOCH, DS50
DATA RHO/.15696615/
IF (IFLAG .EQ. 0) GO TO 100
* RECOVER ORIGINAL MEAN MOTION (XNODP) AND SEMIMAJOR AXIS (AODP)
* FROM INPUT ELEMENTS --------- CALCULATE BALLISTIC COEFFICIENT
* (B TERM) FROM INPUT B* DRAG TERM
A1=(XKE/XNO)**TOTHRD
COSI=COS(XINCL)
THETA2=COSI*COSI
TTHMUN=3.*THETA2-1.
EOSQ=EO*EO
BETAO2=1.-EOSQ
BETAO=SQRT(BETAO2)
DEL1=1.5*CK2*TTHMUN/(A1*A1*BETAO*BETAO2)
AO=A1*(1.-DEL1*(.5*TOTHRD+DEL1*(1.+134./81.*DEL1)))
DELO=1.5*CK2*TTHMUN/(AO*AO*BETAO*BETAO2)
AODP=AO/(1.-DELO)
XNODP=XNO/(1.+DELO)
B=2.*BSTAR/RHO
* INITIALIZATION
ISIMP=0
PO=AODP*BETAO2
POM2=1./(PO*PO)
SINI=SIN(XINCL)
SING=SIN(OMEGAO)
COSG=COS(OMEGAO)
TEMP=.5*XINCL
SINIO2=SIN(TEMP)
COSIO2=COS(TEMP)
THETA4=THETA2**2
UNM5TH=1.-5.*THETA2
UNMTH2=1.-THETA2
A3COF=-XJ3/CK2*AE**3
PARDT1=3.*CK2*POM2*XNODP
PARDT2=PARDT1*CK2*POM2
PARDT4=1.25*CK4*POM2*POM2*XNODP
XMDT1=.5*PARDT1*BETAO*TTHMUN
XGDT1=-.5*PARDT1*UNM5TH
XHDT1=-PARDT1*COSI
XLLDOT=XNODP+XMDT1+
2 .0625*PARDT2*BETAO*(13.-78.*THETA2+137.*THETA4)
OMGDT=XGDT1+
1 .0625*PARDT2*(7.-114.*THETA2+395.*THETA4)+PARDT4*(3.-36.*
2 THETA2+49.*THETA4)
XNODOT=XHDT1+
1 (.5*PARDT2*(4.-19.*THETA2)+2.*PARDT4*(3.-7.*THETA2))*COSI
TSI=1./(PO-S)
ETA=EO*S*TSI
ETA2=ETA**2
PSIM2=ABS(1./(1.-ETA2))
ALPHA2=1.+EOSQ
EETA=EO*ETA
COS2G=2.*COSG**2-1.
D5=TSI*PSIM2
D1=D5/PO
D2=12.+ETA2*(36.+4.5*ETA2)
D3=ETA2*(15.+2.5*ETA2)
D4=ETA*(5.+3.75*ETA2)
B1=CK2*TTHMUN
B2=-CK2*UNMTH2
B3=A3COF*SINI
C0=.5*B*RHO*QOMS2T*XNODP*AODP*TSI**4*PSIM2**3.5/SQRT(ALPHA2)
C1=1.5*XNODP*ALPHA2**2*C0
C4=D1*D3*B2
C5=D5*D4*B3
XNDT=C1*(
1 (2.+ETA2*(3.+34.*EOSQ)+5.*EETA*(4.+ETA2)+8.5*EOSQ)+
1 D1*D2*B1+ C4*COS2G+C5*SING)
XNDTN=XNDT/XNODP
* IF DRAG IS VERY SMALL, THE ISIMP FLAG IS SET AND THE
* EQUATIONS ARE TRUNCATED TO LINEAR VARIATION IN MEAN
* MOTION AND QUADRATIC VARIATION IN MEAN ANOMALY
IF(ABS(XNDTN*XMNPDA) .LT. 2.16E-3) GO TO 50
D6=ETA*(30.+22.5*ETA2)
D7=ETA*(5.+12.5*ETA2)
D8=1.+ETA2*(6.75+ETA2)
C8=D1*D7*B2
C9=D5*D8*B3
EDOT=-C0*(
1 ETA*(4.+ETA2+EOSQ*(15.5+7.*ETA2))+EO*(5.+15.*ETA2)+
1 D1*D6*B1 +
1 C8*COS2G+C9*SING)
D20=.5*TOTHRD*XNDTN
ALDTAL=EO*EDOT/ALPHA2
TSDTTS=2.*AODP*TSI*(D20*BETAO2+EO*EDOT)
ETDT=(EDOT+EO*TSDTTS)*TSI*S
PSDTPS=-ETA*ETDT*PSIM2
SIN2G=2.*SING*COSG
C0DTC0=D20+4.*TSDTTS-ALDTAL-7.*PSDTPS
C1DTC1=XNDTN+4.*ALDTAL+C0DTC0
D9=ETA*(6.+68.*EOSQ)+EO*(20.+15.*ETA2)
D10=5.*ETA*(4.+ETA2)+EO*(17.+68.*ETA2)
D11=ETA*(72.+18.*ETA2)
D12=ETA*(30.+10.*ETA2)
D13=5.+11.25*ETA2
D14=TSDTTS-2.*PSDTPS
D15=2.*(D20+EO*EDOT/BETAO2)
D1DT=D1*(D14+D15)
D2DT=ETDT*D11
D3DT=ETDT*D12
D4DT=ETDT*D13
D5DT=D5*D14
C4DT=B2*(D1DT*D3+D1*D3DT)
C5DT=B3*(D5DT*D4+D5*D4DT)
D16=
1 D9*ETDT+D10*EDOT +
1 B1*(D1DT*D2+D1*D2DT) +
1 C4DT*COS2G+C5DT*SING+XGDT1*(C5*COSG-2.*C4*SIN2G)
XNDDT=C1DTC1*XNDT+C1*D16
EDDOT=C0DTC0*EDOT-C0*(
1 (4.+3.*ETA2+30.*EETA+EOSQ*(15.5+21.*ETA2))*ETDT+(5.+15.*ETA2
' +EETA*(31.+14.*ETA2))*EDOT +
1 B1*(D1DT*D6+D1*ETDT*(30.+67.5*ETA2)) +
1 B2*(D1DT*D7+D1*ETDT*(5.+37.5*ETA2))*COS2G+
1 B3*(D5DT*D8+D5*ETDT*ETA*(13.5+4.*ETA2))*SING+XGDT1*(C9*
' COSG-2.*C8*SIN2G))
D25=EDOT**2
D17=XNDDT/XNODP-XNDTN**2
TSDDTS=2.*TSDTTS*(TSDTTS-D20)+AODP*TSI*(TOTHRD*BETAO2*D17-4.*D20*
' EO*EDOT+2.*(D25+EO*EDDOT))
ETDDT =(EDDOT+2.*EDOT*TSDTTS)*TSI*S+TSDDTS*ETA
D18=TSDDTS-TSDTTS**2
D19=-PSDTPS**2/ETA2-ETA*ETDDT*PSIM2-PSDTPS**2
D23=ETDT*ETDT
D1DDT=D1DT*(D14+D15)+D1*(D18-2.*D19+TOTHRD*D17+2.*(ALPHA2*D25
' /BETAO2+EO*EDDOT)/BETAO2)
XNTRDT=XNDT*(2.*TOTHRD*D17+3.*
1 (D25+EO*EDDOT)/ALPHA2-6.*ALDTAL**2 +
1 4.*D18-7.*D19 ) +
1 C1DTC1*XNDDT+C1*(C1DTC1*D16+
1 D9*ETDDT+D10*EDDOT+D23*(6.+30.*EETA+68.*EOSQ)+
1 ETDT*EDOT*(40.+30.*
' ETA2+272.*EETA)+D25*(17.+68.*ETA2) +
1 B1*(D1DDT*D2+2.*D1DT*D2DT+D1*(ETDDT*D11+D23*(72.+54.*ETA2))) +
1 B2*(D1DDT*D3+2.*D1DT*D3DT+D1*(ETDDT*D12+D23*(30.+30.*ETA2))) *
1 COS2G+
1 B3*((D5DT*D14+D5*(D18-2.*D19)) *
1 D4+2.*D4DT*D5DT+D5*(ETDDT*D13+22.5*ETA*D23)) *SING+XGDT1*
1 ((7.*D20+4.*EO*EDOT/BETAO2)*
' (C5*COSG-2.*C4*SIN2G)
' +((2.*C5DT*COSG-4.*C4DT*SIN2G)-XGDT1*(C5*SING+4.*
' C4*COS2G))))
TMNDDT=XNDDT*1.E9
TEMP=TMNDDT**2-XNDT*1.E18*XNTRDT
PP=(TEMP+TMNDDT**2)/TEMP
GAMMA=-XNTRDT/(XNDDT*(PP-2.))
XND=XNDT/(PP*GAMMA)
QQ=1.-EDDOT/(EDOT*GAMMA)
ED=EDOT/(QQ*GAMMA)
OVGPP=1./(GAMMA*(PP+1.))
GO TO 70
50 ISIMP=1
EDOT=-TOTHRD*XNDTN*(1.-EO)
70 IFLAG=0
* UPDATE FOR SECULAR GRAVITY AND ATMOSPHERIC DRAG
100 XMAM=FMOD2P(XMO+XLLDOT*TSINCE)
OMGASM=OMEGAO+OMGDT*TSINCE
XNODES=XNODEO+XNODOT*TSINCE
IF(ISIMP .EQ. 1) GO TO 105
TEMP=1.-GAMMA*TSINCE
TEMP1=TEMP**PP
XN=XNODP+XND*(1.-TEMP1)
EM=EO+ED*(1.-TEMP**QQ)
Z1=XND*(TSINCE+OVGPP*(TEMP*TEMP1-1.))
GO TO 108
105 XN=XNODP+XNDT*TSINCE
EM=EO+EDOT*TSINCE
Z1=.5*XNDT*TSINCE*TSINCE
108 Z7=3.5*TOTHRD*Z1/XNODP
XMAM=FMOD2P(XMAM+Z1+Z7*XMDT1)
OMGASM=OMGASM+Z7*XGDT1
XNODES=XNODES+Z7*XHDT1
* SOLVE KEPLERS EQUATION
ZC2=XMAM+EM*SIN(XMAM)*(1.+EM*COS(XMAM))
DO 130 I=1,10
SINE=SIN(ZC2)
COSE=COS(ZC2)
ZC5=1./(1.-EM*COSE)
CAPE=(XMAM+EM*SINE-ZC2)*
1 ZC5+ZC2
IF(ABS(CAPE-ZC2) .LE. E6A) GO TO 140
130 ZC2=CAPE
* SHORT PERIOD PRELIMINARY QUANTITIES
140 AM=(XKE/XN)**TOTHRD
BETA2M=1.-EM*EM
SINOS=SIN(OMGASM)
COSOS=COS(OMGASM)
AXNM=EM*COSOS
AYNM=EM*SINOS
PM=AM*BETA2M
G1=1./PM
G2=.5*CK2*G1
G3=G2*G1
BETA=SQRT(BETA2M)
G4=.25*A3COF*SINI
G5=.25*A3COF*G1
SNF=BETA*SINE*ZC5
CSF=(COSE-EM)*ZC5
FM=ACTAN(SNF,CSF)
SNFG=SNF*COSOS+CSF*SINOS
CSFG=CSF*COSOS-SNF*SINOS
SN2F2G=2.*SNFG*CSFG
CS2F2G=2.*CSFG**2-1.
ECOSF=EM*CSF
G10=FM-XMAM+EM*SNF
RM=PM/(1.+ECOSF)
AOVR=AM/RM
G13=XN*AOVR
G14=-G13*AOVR
DR=G2*(UNMTH2*CS2F2G-3.*TTHMUN)-G4*SNFG
DIWC=3.*G3*SINI*CS2F2G-G5*AYNM
DI=DIWC*COSI
* UPDATE FOR SHORT PERIOD PERIODICS
SNI2DU=SINIO2*(
1 G3*(.5*(1.-7.*THETA2)*SN2F2G-3.*UNM5TH*G10)-G5*SINI*CSFG*(2.+
2 ECOSF))-.5*G5*THETA2*AXNM/COSIO2
XLAMB=FM+OMGASM+XNODES+G3*(.5*(1.+6.*COSI-7.*THETA2)*SN2F2G-3.*
1 (UNM5TH+2.*COSI)*G10)+G5*SINI*(COSI*AXNM/(1.+COSI)-(2.
2 +ECOSF)*CSFG)
Y4=SINIO2*SNFG+CSFG*SNI2DU+.5*SNFG*COSIO2*DI
Y5=SINIO2*CSFG-SNFG*SNI2DU+.5*CSFG*COSIO2*DI
R=RM+DR
RDOT=XN*AM*EM*SNF/BETA+G14*(2.*G2*UNMTH2*SN2F2G+G4*CSFG)
RVDOT=XN*AM**2*BETA/RM+
1 G14*DR+AM*G13*SINI*DIWC
* ORIENTATION VECTORS
SNLAMB=SIN(XLAMB)
CSLAMB=COS(XLAMB)
TEMP=2.*(Y5*SNLAMB-Y4*CSLAMB)
UX=Y4*TEMP+CSLAMB
VX=Y5*TEMP-SNLAMB
TEMP=2.*(Y5*CSLAMB+Y4*SNLAMB)
UY=-Y4*TEMP+SNLAMB
VY=-Y5*TEMP+CSLAMB
TEMP=2.*SQRT(1.-Y4*Y4-Y5*Y5)
UZ=Y4*TEMP
VZ=Y5*TEMP
* POSITION AND VELOCITY
X=R*UX
Y=R*UY
Z=R*UZ
XDOT=RDOT*UX+RVDOT*VX
YDOT=RDOT*UY+RVDOT*VY
ZDOT=RDOT*UZ+RVDOT*VZ
RETURN
END
– SDP8.FOR –
* SDP8 14 NOV 80
SUBROUTINE SDP8(IFLAG,TSINCE)
COMMON/E1/XMO,XNODEO,OMEGAO,EO,XINCL,XNO,XNDT2O,
1 XNDD6O,BSTAR,X,Y,Z,XDOT,YDOT,ZDOT,EPOCH,DS50
COMMON/C1/CK2,CK4,E6A,QOMS2T,S,TOTHRD,
1 XJ3,XKE,XKMPER,XMNPDA,AE
DOUBLE PRECISION EPOCH, DS50
DATA RHO/.15696615/
IF (IFLAG .EQ. 0) GO TO 100
* RECOVER ORIGINAL MEAN MOTION (XNODP) AND SEMIMAJOR AXIS (AODP)
* FROM INPUT ELEMENTS --------- CALCULATE BALLISTIC COEFFICIENT
* (B TERM) FROM INPUT B* DRAG TERM
A1=(XKE/XNO)**TOTHRD
COSI=COS(XINCL)
THETA2=COSI*COSI
TTHMUN=3.*THETA2-1.
EOSQ=EO*EO
BETAO2=1.-EOSQ
BETAO=SQRT(BETAO2)
DEL1=1.5*CK2*TTHMUN/(A1*A1*BETAO*BETAO2)
AO=A1*(1.-DEL1*(.5*TOTHRD+DEL1*(1.+134./81.*DEL1)))
DELO=1.5*CK2*TTHMUN/(AO*AO*BETAO*BETAO2)
AODP=AO/(1.-DELO)
XNODP=XNO/(1.+DELO)
B=2.*BSTAR/RHO
* INITIALIZATION
PO=AODP*BETAO2
POM2=1./(PO*PO)
SINI=SIN(XINCL)
SING=SIN(OMEGAO)
COSG=COS(OMEGAO)
TEMP=.5*XINCL
SINIO2=SIN(TEMP)
COSIO2=COS(TEMP)
THETA4=THETA2**2
UNM5TH=1.-5.*THETA2
UNMTH2=1.-THETA2
A3COF=-XJ3/CK2*AE**3
PARDT1=3.*CK2*POM2*XNODP
PARDT2=PARDT1*CK2*POM2
PARDT4=1.25*CK4*POM2*POM2*XNODP
XMDT1=.5*PARDT1*BETAO*TTHMUN
XGDT1=-.5*PARDT1*UNM5TH
XHDT1=-PARDT1*COSI
XLLDOT=XNODP+XMDT1+
2 .0625*PARDT2*BETAO*(13.-78.*THETA2+137.*THETA4)
OMGDT=XGDT1+
1 .0625*PARDT2*(7.-114.*THETA2+395.*THETA4)+PARDT4*(3.-36.*
2 THETA2+49.*THETA4)
XNODOT=XHDT1+
1 (.5*PARDT2*(4.-19.*THETA2)+2.*PARDT4*(3.-7.*THETA2))*COSI
TSI=1./(PO-S)
ETA=EO*S*TSI
ETA2=ETA**2
PSIM2=ABS(1./(1.-ETA2))
ALPHA2=1.+EOSQ
EETA=EO*ETA
COS2G=2.*COSG**2-1.
D5=TSI*PSIM2
D1=D5/PO
D2=12.+ETA2*(36.+4.5*ETA2)
D3=ETA2*(15.+2.5*ETA2)
D4=ETA*(5.+3.75*ETA2)
B1=CK2*TTHMUN
B2=-CK2*UNMTH2
B3=A3COF*SINI
C0=.5*B*RHO*QOMS2T*XNODP*AODP*TSI**4*PSIM2**3.5/SQRT(ALPHA2)
C1=1.5*XNODP*ALPHA2**2*C0
C4=D1*D3*B2
C5=D5*D4*B3
XNDT=C1*(
1 (2.+ETA2*(3.+34.*EOSQ)+5.*EETA*(4.+ETA2)+8.5*EOSQ)+
1 D1*D2*B1+ C4*COS2G+C5*SING)
XNDTN=XNDT/XNODP
EDOT=-TOTHRD*XNDTN*(1.-EO)
IFLAG=0
CALL DPINIT(EOSQ,SINI,COSI,BETAO,AODP,THETA2,SING,COSG,
1 BETAO2,XLLDOT,OMGDT,XNODOT,XNODP)
* UPDATE FOR SECULAR GRAVITY AND ATMOSPHERIC DRAG
100 Z1=.5*XNDT*TSINCE*TSINCE
Z7=3.5*TOTHRD*Z1/XNODP
XMAMDF=XMO+XLLDOT*TSINCE
OMGASM=OMEGAO+OMGDT*TSINCE+Z7*XGDT1
XNODES=XNODEO+XNODOT*TSINCE+Z7*XHDT1
XN=XNODP
CALL DPSEC(XMAMDF,OMGASM,XNODES,EM,XINC,XN,TSINCE)
XN=XN+XNDT*TSINCE
EM=EM+EDOT*TSINCE
XMAM=XMAMDF+Z1+Z7*XMDT1
CALL DPPER(EM,XINC,OMGASM,XNODES,XMAM)
XMAM=FMOD2P(XMAM)
* SOLVE KEPLERS EQUATION
ZC2=XMAM+EM*SIN(XMAM)*(1.+EM*COS(XMAM))
DO 130 I=1,10
SINE=SIN(ZC2)
COSE=COS(ZC2)
ZC5=1./(1.-EM*COSE)
CAPE=(XMAM+EM*SINE-ZC2)*
1 ZC5+ZC2
IF(ABS(CAPE-ZC2) .LE. E6A) GO TO 140
130 ZC2=CAPE
* SHORT PERIOD PRELIMINARY QUANTITIES
140 AM=(XKE/XN)**TOTHRD
BETA2M=1.-EM*EM
SINOS=SIN(OMGASM)
COSOS=COS(OMGASM)
AXNM=EM*COSOS
AYNM=EM*SINOS
PM=AM*BETA2M
G1=1./PM
G2=.5*CK2*G1
G3=G2*G1
BETA=SQRT(BETA2M)
G4=.25*A3COF*SINI
G5=.25*A3COF*G1
SNF=BETA*SINE*ZC5
CSF=(COSE-EM)*ZC5
FM=ACTAN(SNF,CSF)
SNFG=SNF*COSOS+CSF*SINOS
CSFG=CSF*COSOS-SNF*SINOS
SN2F2G=2.*SNFG*CSFG
CS2F2G=2.*CSFG**2-1.
ECOSF=EM*CSF
G10=FM-XMAM+EM*SNF
RM=PM/(1.+ECOSF)
AOVR=AM/RM
G13=XN*AOVR
G14=-G13*AOVR
DR=G2*(UNMTH2*CS2F2G-3.*TTHMUN)-G4*SNFG
DIWC=3.*G3*SINI*CS2F2G-G5*AYNM
DI=DIWC*COSI
SINI2=SIN(.5*XINC)
* UPDATE FOR SHORT PERIOD PERIODICS
SNI2DU=SINIO2*(
1 G3*(.5*(1.-7.*THETA2)*SN2F2G-3.*UNM5TH*G10)-G5*SINI*CSFG*(2.+
2 ECOSF))-.5*G5*THETA2*AXNM/COSIO2
XLAMB=FM+OMGASM+XNODES+G3*(.5*(1.+6.*COSI-7.*THETA2)*SN2F2G-3.*
1 (UNM5TH+2.*COSI)*G10)+G5*SINI*(COSI*AXNM/(1.+COSI)-(2.
2 +ECOSF)*CSFG)
Y4=SINI2*SNFG+CSFG*SNI2DU+.5*SNFG*COSIO2*DI
Y5=SINI2*CSFG-SNFG*SNI2DU+.5*CSFG*COSIO2*DI
R=RM+DR
RDOT=XN*AM*EM*SNF/BETA+G14*(2.*G2*UNMTH2*SN2F2G+G4*CSFG)
RVDOT=XN*AM**2*BETA/RM+
1 G14*DR+AM*G13*SINI*DIWC
* ORIENTATION VECTORS
SNLAMB=SIN(XLAMB)
CSLAMB=COS(XLAMB)
TEMP=2.*(Y5*SNLAMB-Y4*CSLAMB)
UX=Y4*TEMP+CSLAMB
VX=Y5*TEMP-SNLAMB
TEMP=2.*(Y5*CSLAMB+Y4*SNLAMB)
UY=-Y4*TEMP+SNLAMB
VY=-Y5*TEMP+CSLAMB
TEMP=2.*SQRT(1.-Y4*Y4-Y5*Y5)
UZ=Y4*TEMP
VZ=Y5*TEMP
* POSITION AND VELOCITY
X=R*UX
Y=R*UY
Z=R*UZ
XDOT=RDOT*UX+RVDOT*VX
YDOT=RDOT*UY+RVDOT*VY
ZDOT=RDOT*UZ+RVDOT*VZ
RETURN
END
– THETAG.FOR –
FUNCTION THETAG(EP)
COMMON /E1/XMO,XNODEO,OMEGAO,EO,XINCL,XNO,XNDT2O,XNDD6O,BSTAR,
1 X,Y,Z,XDOT,YDOT,ZDOT,EPOCH,DS50
DOUBLE PRECISION EPOCH,D,THETA,TWOPI,YR,TEMP,EP,DS50
TWOPI=6.28318530717959D0
YR=(EP+2.D-7)*1.D-3
JY=YR
YR=JY
D=EP-YR*1.D3
IF(JY.LT.10) JY=JY+80
N=(JY-69)/4
IF(JY.LT.70) N=(JY-72)/4
DS50=7305.D0 + 365.D0*(JY-70) +N + D
THETA=1.72944494D0 + 6.3003880987D0*DS50
TEMP=THETA/TWOPI
I=TEMP
TEMP=I
THETAG=THETA-TEMP*TWOPI
IF(THETAG.LT.0.D0) THETAG=THETAG+TWOPI
RETURN
END